#
rm(list=ls())
library(tictoc)
library(HeadR)
library(doParallel)
library(SQUAREM)
library(fixest)
setwd("~/Dropbox/BLP_REStat/Replication")
library(Rcpp)
# create the Rcpp function we need for computing equilibrium
Rcpp::sourceCpp("Cppfiles/crossDelta.cpp")
#parallel processing set up for doParallel, doRNG
cl <- detectCores()-1
registerDoParallel(cores=cl)
getDoParWorkers()
MAXIT <- 500
seedn <- 777 # 666
sequences <- "normal" # alternatives: normal, halton or "sobol"
sigmadomesticpct <- 0
#
n.hh <- 100000 # careful, this is slow, but works
U0 <- 1  # needed in prob.mat(), but not prob.mat.namedxvars()
n.reps <- 1 # for richsubs (no need to do a lot if a lot of consumers.)
speed.list <- seq(0.2,1,by=0.05) # different dampenings
num.loops <- length(speed.list)
#
# This is intended to choose which country is targeted
source(file = "Rfiles/CF_BLPdata_brandiso.R")
target.list <- brands[iso!="USA",unique(iso)] # all non US brands
rm(brands)
#
#
# Here choose which set of parameters to use. 
param.all <- fread("Data/BLP99/pyblp_domestic_param.csv") # estimated using C&G package
param.list <- names(param.all[,par_abbrev:=NULL])
param.list 
#add pub_param (the BLP 99 parameters)
param.chosen <- c("pub_param")
#
#
tic()
##
#Loop over sigmadomesticpct list, conducting counterfactuals ====
##
DF <- foreach(k=1:num.loops, .combine = 'rbind') %dopar% {
#
source("Rfiles/CF_BLPdata_domestic_functions.R",local=TRUE)
source("Rfiles/CF_MMX_functions.R",local=TRUE)
source("Rfiles/CF_MMX_MLOG_functions.R",local=TRUE)
#  
settings.var <- data.frame(speed =rep(speed.list[k],3))
#
cat("Now running for speed  =", speed.list[k] ,"\n")
#
# This is used in the counterfactuals : number of characteristics (outside the constant)
if(param.chosen %like% "^Domestic"){
n.x <- 5 #THIS NEEDS TO BE 5 WHEN ADDING DOMESTIC, 4 otherwise
} else n.x <- 4
#
###
# What is the average own elasticity? ====
###
setwd("~/Dropbox/BLP_REStat/Replication")
source("Rfiles/CF_BLPdata_domestic_getdata.R",local = TRUE) # get original data: original published_param.csv is 1995 and should not be changed (used in other files)
#
year.gph <- 1990
year.monte <- year.gph
ybar <- exp(y.mn[year ==year.gph]$V2)
# co ownership matrix stuff
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
#initialization phase
set.seed(seedn)
data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=3) #gen exog data with v=3 for initialization
BLP.output <- rich.subs.BLP.mlog(dat=data.exog.mm,verbose=FALSE,objective="richsubs") 
BLP.output.m <- BLP.output[m==j]
BLP.meanownelas <- BLP.output.m[, mean(own_elas_m)] # super assignment needed here because of scope for DGP.BLP
BLP.meanownelas
#
###
# Implement CF on BLP ====
###
#
source("Rfiles/CF_BLPdata_domestic_getdata.R",local = TRUE) # get original data
year.gph <- 1990
tar <- 0 # original tariff (to be raised)
tar.increase <- 0.1 # the actual increase in tariff rate
ybar <- exp(y.mn[year ==year.gph]$V2)
# co ownership matrix 
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
#
# Run the calibration to get xi and mc
year.monte <- year.gph
objective <- "CF"
market <- 1
#
#Run sims for the 2 relevant versions (no logit)
#EHA:
set.seed(seedn)
CF_BLPdata.simulv3 <- CF_BLPdata.uniqueness(vsim=3,cf.meth="EHA")
CF_BLPdata.simulv3[, speed := speed.list[k]]
#
CF_BLPdata.simulv3
}
toc()
#
#
DFU <- DF[, .(s_sd = sd(s), max.loops = max(n.loops)) , by=m]
DFU[, .(max_sd=max(s_sd),max.nloops=max(max.loops))]
DF[, .(max_Delta=max(chg.true.qty),min_Delta=min(chg.true.qty))]
#
plot(unique(DF[, .(speed,n.loops)]),type="b",pch=15)
